c2.rat=function(data,ybari,s2i,mi,Mi,n,N,param){
  data=as.name(data); N=N; mi=mi; n=n
  Mi=Mi; fpc=1-n/N; ybari=ybari
  s2i=s2i; Mbar=sum(Mi)/n
  muhatr=sum(Mi*ybari)/sum(Mi)
  s2r=sum((Mi*ybari-Mi*muhatr)^2)/(n-1)
  if(param=='mean'){
    est=muhatr
    vhat=fpc*s2r*(1/(n*Mbar^2))+(1/(n*N*Mbar^2))*sum(Mi^2*(1-mi/Mi)*(s2i/mi))
    bound=2*sqrt(vhat)
    lower=est-bound; upper=est+bound
  } else if(param=='total'){
    est=N*muhatr
    vhat=N^2*(fpc*s2r*(1/(n*Mbar^2))+(1/(n*N*Mbar^2))*sum(Mi^2*(1-mi/Mi)*(s2i/mi)))
    bound=2*sqrt(vhat)
    lower=est-bound; upper=est+bound
  } 
  cat("","\n","Results from 2-stage Cluster design: Data =",data,"\n","\n","N =",N,"n =",n,"FPC =",fpc,"\n",
      "Estimate of",param,"=",est,"\n","Vhat(", param,") =",vhat,"\n","Bound =",bound,"\n","Lower Bound =",
      lower,"Upper Bound =",upper,"\n","")
  results=list(est=est,data=data,n=n,N=N,fpc=fpc,vhat=vhat,bound=bound,lower=lower,upper=upper,param=param,
               s2i=s2i,s2r=s2r)
}

# to use the function with its call:
# c2.rat(data,ybari,s2i,mi,Mi,n,N,param)

# data: name of dataset, in quotes
# ybari: numeric vector of the ith cluster means
# s2i: numeric vector of the ith cluster variances
# mi: numeric vector of the ith cluster elements
# Mi: numeric vector of the population cluster elements
# N: single number of population size
# param: c('mean','total')